ftp.mactech.com 2010/ftp.mactech.com.tar/ftp.mactech.com/challenge/13.07/Challenge.sit/Challenge, Polygon Projection/Challenge.p

Text File | 1997-05-07 | 16KB | 485 lines

unit Challenge;

(*

    Author: A.C.C. Murphy

    Assumptions:
        Storage space must be big enough for 13 floats per polygon
        All points must be significantly smaller in magnitude than BIG_FLOAT = 1000000.0
        Polygons are translucent (their colour based uplon lighting is independent of the side of the
            polygon that is lit)
        50% attenuation of colour is used
        50% attenuation of black is black

    Method:

        InitProjection is not used

        First we precalculate a small bounding sphere for the polygon points.
        Next we get the information about the GWorld to allow direct pixel access.
        Then for each point on the GWorld, we trace the ray from the point to the eye,
            intersecting it with each polygon and finding the one closes to the eye
            (furthest forward, since the eye is infront of all polygons).  That determines
            the colour.  We then trace the ray from that intersection point to the light source
            to determine whether the point is in shadow, and if so we halve the intensity.
            We set the colour of the pixel and move on.

    Optimizations:
        Direct pixel access to the GWorld (known to be 32 bit)
        Bounding sphere used to optimize the ray/polygon intersection test.
        Time is approximately 2 microseconds per pixel per polygon on an 8500.
*)

interface

    uses
        Types, Quickdraw, QDOffscreen;

    const
        kMAXPOINTS = 10;

    const
        BIG_FLOAT = 1000000.0;

    type
        float = real;

    type
        My2DPoint = record (* point in z==0 plane*)
            x2D: float; (* x coordinate*)
            y2D: float; (* y coordinate*)
        end;
        My3DPoint = record
            x3D: float;                 (* x coordinate*)
            y3D: float;                 (* y coordinate*)
            z3D: float;                 (* z coordinate*)
        end;
        My3DDirection = record
            thetaX:float;              (* angle in radians*)
            thetaY:float;              (* angle in radians*)
            thetaZ:float;              (* angle in radians*)
        end;
        MyPlane = record
            planeNormal: My3DDirection;      (* normal vector to plane*)
            planeOrigin: My3DPoint;    (* origin of plane in 3D space*)
        end;
        MyPolygon = record
            numPoints: longint;             (* number of points in polygon*)
            thePoint: array[0..kMAXPOINTS-1] of My2DPoint;     (* polygon in z==0 plane*)
            polyPlane: MyPlane;             (* rotate/translate z==0 plane to this plane*)
            polyColor: RGBColor;         (* the color to draw this polygon*)
        end;
        MyPolygonArray = array[0..0] of MyPolygon;


    procedure InitProjection(
        const viewPoint: My3DPoint;         (* viewpoint from which to project*)
        const illumPoint:My3DPoint;         (* viewpoint from which to draw shadow*)
        storage: univ Ptr;                         (* auxiliary storage preallocated for your use*)
        storageSize: longint                     (* number of bytes of storage*)
    );

    procedure CalcProjection(
        offScreen: GWorldPtr;                                 (* GWorld to draw projection *)
        const thePolys: MyPolygonArray;         (* polygons to project *)
        numPolys: longint;                                         (* number of polygons to project *)
        const viewPoint: My3DPoint;                             (* viewpoint from which to project *)
        const illumPoint: My3DPoint;                             (* illumination point from which to draw shadow *)
        storage: univ Ptr;                         (* auxiliary storage preallocated for your use*)
        storageSize: longint                     (* number of bytes of storage*)
    );

implementation

    type
        Ray3D = record
            origin: My3DPoint;
            direction: My3DPoint;
        end;
        PolygonExtra = record
            normal, rotX, rotY, center: My3DPoint;
            radius2: float;
        end;
        PolygonExtraArray = array[0..0] of PolygonExtra;
        StorageRecord = record
            poly_extra: PolygonExtraArray; { must be at the end, since it's an extensible array }
        end;
        StorageRecordPtr = ^StorageRecord;

    function DotProduct(const src1, src2 : My3DPoint) : float;
    begin
        DotProduct := src1.x3D*src2.x3D +  src1.y3D*src2.y3D +  src1.z3D*src2.z3D;
    end;

    procedure CrossProduct(src1, src2 : My3DPoint; var dst : My3DPoint);
    begin
        dst.x3D := src1.y3D*src2.z3D - src1.z3D*src2.y3D;
        dst.y3D := src1.z3D*src2.x3D - src1.x3D*src2.z3D;
        dst.z3D := src1.x3D*src2.y3D - src1.y3D*src2.x3D;
    end;

    procedure AddVectors(const src1, src2 : My3DPoint; var dst : My3DPoint);
    begin
        dst.x3D := src1.x3D + src2.x3D;
        dst.y3D := src1.y3D + src2.y3D;
        dst.z3D := src1.z3D + src2.z3D;
    end;

    procedure SubtractVectors(const src1, src2 : My3DPoint; var dst : My3DPoint);
    begin
        dst.x3D := src1.x3D - src2.x3D;
        dst.y3D := src1.y3D - src2.y3D;
        dst.z3D := src1.z3D - src2.z3D;
    end;

    procedure MidPoint( const src1, src2 : My3DPoint; var dst : My3DPoint);
    begin
        dst.x3D := (src1.x3D + src2.x3D) / 2;
        dst.y3D := (src1.y3D + src2.y3D) / 2;
        dst.z3D := (src1.z3D + src2.z3D) / 2;
    end;

    function Distance2( const src1, src2 : My3DPoint) : float;
    begin
        Distance2 := sqr(src1.x3D - src2.x3D) + sqr(src1.y3D - src2.y3D) + sqr(src1.z3D - src2.z3D);
    end;

    procedure ScaleVector(const src : My3DPoint; scale : float; var dst : My3DPoint);
    begin
        dst.x3D := src.x3D * scale;
        dst.y3D := src.y3D * scale;
        dst.z3D := src.z3D * scale;
    end;

    procedure NormalizeVector(const src : My3DPoint; var dst : My3DPoint);
        var
            length : float;
    begin
        length := sqrt(DotProduct(src,src));
        dst.x3D := src.x3D / length;
        dst.y3D := src.y3D / length;
        dst.z3D := src.z3D / length;
    end;

    procedure MakeViewRay(const eye : My3DPoint; x, y, z: float; var ray : Ray3D);
    begin
        ray.origin.x3D := x;
        ray.origin.y3D := y;
        ray.origin.z3D := z;
        ray.direction.x3D := eye.x3D - x;
        ray.direction.y3D := eye.y3D - y;
        ray.direction.z3D := eye.z3D - z;
        NormalizeVector(ray.direction, ray.direction);
    end;

    procedure RotateX(src : My3DPoint; sinA, cosA : float; var dst : My3DPoint);
    begin
        dst.x3D := src.x3D;
        dst.y3D := cosA*src.y3D - sinA*src.z3D;
        dst.z3D := sinA*src.y3D + cosA*src.z3D;
    end;

    procedure RotateY( src : My3DPoint; sinA, cosA : float; var dst : My3DPoint);
    begin
        dst.x3D := cosA*src.x3D + sinA*src.z3D;
        dst.y3D := src.y3D;
        dst.z3D := -sinA*src.x3D + cosA*src.z3D;
    end;

    procedure RotateZ( src : My3DPoint; sinA, cosA : float; var dst : My3DPoint);
    begin
        dst.x3D := cosA*src.x3D - sinA*src.y3D;
        dst.y3D := sinA*src.x3D + cosA*src.y3D;
        dst.z3D := src.z3D;
    end;

    function PointInPlaneInPolygon( const pt: My2DPoint; const poly: MyPolygon ): boolean;
        function Quadrant( const pt: My2DPoint; x, y: float ): longint;
        begin
            if pt.x2D > x then begin
                if pt.y2D > y then begin
                    Quadrant := 0;
                end else begin
                    Quadrant := 3;
                end;
            end else begin
                if pt.y2D > y then begin
                    Quadrant := 1;
                end else begin
                    Quadrant := 2;
                end;
            end;
        end;

        function x_intercept( const pt1, pt2: My2DPoint; yy: float ): float;
        begin
            x_intercept := pt2.x2D - ( (pt2.y2D - yy) * ((pt1.x2D - pt2.x2D) / (pt1.y2D - pt2.y2D)) );
        end;

        var
            i, angle, quad, next_quad, delta: longint;
            last_vertex, next_vertex: My2DPoint;
    begin
        angle := 0;
        last_vertex := poly.thePoint[poly.numPoints-1];
        quad := Quadrant( last_vertex, pt.x2D, pt.y2D );
        for i := 1 to poly.numPoints do begin
            next_vertex := poly.thePoint[i-1];
            next_quad := Quadrant( next_vertex, pt.x2D, pt.y2D );
            delta := next_quad - quad;
            case delta of
                3: delta := -1;
                -3: delta := 1;
                2, -2: begin
                    if x_intercept( last_vertex, next_vertex, pt.y2D ) > pt.x2D then begin
                        delta := -delta;
                    end;
                end;
                otherwise begin
                end;
            end;
            angle := angle + delta;
            quad := next_quad;
            last_vertex := next_vertex;
        end;
        PointInPlaneInPolygon := (angle = 4) | (angle = -4);
    end;

    function Intersect(const ray: Ray3D; const poly: MyPolygon; const poly_extra: PolygonExtra; var distance : float; var intersectionPt: My3DPoint) : boolean;
        var
            tempVector : My3DPoint;
            temp1, temp2 : float;
            intersectionPoint : My3DPoint;
            intersection2D : My2DPoint;
            Ib, Ic, Id: float;
    begin
        Intersect := false;

        { intersect ray with sphere }
        SubtractVectors( ray.origin, poly_extra.center, tempVector );
        Ib := 2 * DotProduct( ray.direction, tempVector );
        Ic := DotProduct( tempVector, tempVector ) - poly_extra.radius2;
        Id := sqr(Ib) - 4.0*Ic;
        if Id >= 0 then begin { yes, ray intersects sphere }
            temp1 := DotProduct( poly.polyPlane.planeOrigin, poly_extra.normal ) - DotProduct( poly_extra.normal, ray.origin );
            temp2 := DotProduct( ray.direction, poly_extra.normal );
            if temp2 <> 0 then begin
                distance := temp1 / temp2;
                if distance > 0 then begin
                    ScaleVector(ray.direction, distance, intersectionPoint);
                    AddVectors(intersectionPoint, ray.origin, intersectionPoint);

                    if Distance2( intersectionPoint, poly_extra.center ) <= poly_extra.radius2 then begin
                        { intersection point is whithin sphere.  Find out if it is actually in the polygon }
                        intersectionPt := intersectionPoint;
                        { First translate back to the origin }
                        SubtractVectors(intersectionPoint, poly.polyPlane.planeOrigin, intersectionPoint);
                        intersection2D.x2D := DotProduct( intersectionPoint, poly_extra.rotX );
                        intersection2D.y2D := DotProduct( intersectionPoint, poly_extra.rotY );
                        { Then check if it is whithin the polygon }
                        Intersect := PointInPlaneInPolygon(intersection2D,poly);
                    end;
                end;
            end;
        end;
    end;

    procedure InitProjection(
        const viewPoint: My3DPoint;        (* viewpoint from which to project *)
        const illumPoint:My3DPoint;        (* viewpoint from which to draw shadow *)
        storage: univ Ptr;                        (* auxiliary storage preallocated for your use *)
        storageSize: longint                    (* number of bytes of storage *)
    );
    begin
{$unused( viewPoint, illumPoint, storage, storageSize )}
    end;

    procedure PreparsePolygons( my_storage: StorageRecordPtr; numPolys: longint; const thePolys: MyPolygonArray );
        var
            i, j: longint;
            pt: My3DPoint;
            pts: array[1..kMAXPOINTS] of My3DPoint;
            min_x, min_y, min_z, max_x, max_y, max_z: My3DPoint;
            dist_x, dist_y, dist_z, new_radius2: float;
            radius, new_radius, old_to_new: float;
            sinX, cosX, sinY, cosY, sinZ, cosZ: float;
    begin
        for i := 0 to numPolys-1 do begin
            with my_storage^.poly_extra[i], thePolys[i], polyPlane.planeNormal do begin
                sinX := sin(thetaX);
                cosX := cos(thetaX);
                sinY := sin(thetaY);
                cosY := cos(thetaY);
                sinZ := sin(thetaZ);
                cosZ := cos(thetaZ);
                normal.x3d := sinY*cosX;
                normal.y3d := -sinX;
                normal.z3d := cosY*cosX;
                rotX.x3D := 1;
                rotX.y3D := 0;
                rotX.z3D := 0;
                RotateZ(rotX, sinZ, cosZ, rotX);
                RotateX(rotX, sinX, cosX, rotX);
                RotateY(rotX, sinY, cosY, rotX);
                rotY.x3D := 0;
                rotY.y3D := 1;
                rotY.z3D := 0;
                RotateZ(rotY, sinZ, cosZ, rotY);
                RotateX(rotY, sinX, cosX, rotY);
                RotateY(rotY, sinY, cosY, rotY);

                for j := 1 to numPoints do begin
                    pt.x3D := thePoint[j-1].x2D;
                    pt.y3D := thePoint[j-1].y2D;
                    pt.z3D := 0;
                    RotateZ(pt, sinZ, cosZ, pt);
                    RotateX(pt, sinX, cosX, pt);
                    RotateY(pt, sinY, cosY, pt);
                    pts[j] := pt;
                    if j = 1 then begin
                        min_x := pt; min_y := pt; min_z := pt;
                        max_x := pt; max_y := pt; max_z := pt;
                    end else begin
                        if pt.x3D < min_x.x3D then begin
                            min_x := pt;
                        end;
                        if pt.y3D < min_y.y3D then begin
                            min_y := pt;
                        end;
                        if pt.z3D < min_z.z3D then begin
                            min_z := pt;
                        end;
                        if pt.x3D > max_x.x3D then begin
                            max_x := pt;
                        end;
                        if pt.y3D > max_y.y3D then begin
                            max_y := pt;
                        end;
                        if pt.z3D > max_z.z3D then begin
                            max_z := pt;
                        end;
                    end;
                end;

                dist_x := Distance2( min_x, max_x );
                dist_y := Distance2( min_y, max_y );
                dist_z := Distance2( min_z, max_z );
                if dist_x > dist_y then begin
                    if dist_x > dist_z then begin
                        radius2 := dist_x/4;
                        MidPoint( min_x, max_x, center );
                    end else begin
                        radius2 := dist_z/4;
                        MidPoint( min_z, max_z, center );
                    end;
                end else begin
                    if dist_y > dist_z then begin
                        radius2 := dist_y/4;
                        MidPoint( min_y, max_y, center );
                    end else begin
                        radius2 := dist_z/4;
                        MidPoint( min_z, max_z, center );
                    end;
                end;

                for j := 1 to numPoints do begin
                    new_radius2 := Distance2( center, pts[j] );
                    if new_radius2 > radius2 then begin
                        radius := sqrt(radius2);
                        new_radius := sqrt(new_radius2);
                        radius2 := (radius + new_radius)/2;
                        old_to_new := radius2 - radius;
                        center.x3D := (radius2*center.x3D + old_to_new*pts[j].x3D)/radius;
                        center.y3D := (radius2*center.y3D + old_to_new*pts[j].y3D)/radius;
                        center.z3D := (radius2*center.z3D + old_to_new*pts[j].z3D)/radius;
                        radius2 := sqr(radius2);
                    end;
                end;

                AddVectors( polyPlane.planeOrigin, center, center );

            end;
        end;
    end;

    procedure CalcProjection(
        offScreen: GWorldPtr;                                (* GWorld to draw projection *)
        const thePolys: MyPolygonArray;        (* polygons to project *)
        numPolys: longint;                                        (* number of polygons to project *)
        const viewPoint: My3DPoint;                            (* viewpoint from which to project *)
        const illumPoint: My3DPoint;                            (* illumination point from which to draw shadow *)
        storage: univ Ptr;                        (* auxiliary storage preallocated for your use *)
        storageSize: longint                    (* number of bytes of storage *)
    );
        var
            bounds: Rect;
            x, y : integer;
            colour : RGBColor;
            viewRay : Ray3D;
            lightRay : Ray3D;
            i : integer;
            closestDistance : float;
            closestIntersectionPt: My3DPoint;
            thisDistance : float;
            intersectionPt: My3DPoint;
            intersect_polygon: longint;
            pm: PixMapHandle;
            junk_boolean: boolean;
            pixels: Ptr;
            rowbytes_add: longint;
            my_storage: StorageRecordPtr;
    begin
{$unused( storage, storageSize )}
        my_storage := StorageRecordPtr(storage);

        PreparsePolygons( my_storage, numPolys, thePolys );

        SetGWorld( offScreen, nil );
        bounds := offScreen^.PortRect;
        pm := GetGWorldPixMap( offScreen );
        junk_boolean := LockPixels( pm );
        pixels := GetPixBaseAddr( pm );
        rowbytes_add := band( pm^^.rowBytes, $3FFF ) - 4 * (bounds.right - bounds.left);

        for y := bounds.top to bounds.bottom-1 do begin
            for x := bounds.left to bounds.right-1 do begin
                MakeViewRay(viewPoint, x, y, 0, viewRay);
                closestDistance := 0.0;
                intersect_polygon := -1;
                for i:= 1 to numPolys do begin
                    if Intersect(viewRay, thePolys[i-1], my_storage^.poly_extra[i-1], thisDistance, intersectionPt) then begin
                        if (thisDistance > closestDistance) then begin
                            intersect_polygon := i;
                            closestDistance := thisDistance;
                            closestIntersectionPt := intersectionPt;
                        end;
                    end
                end;
                if intersect_polygon > 0 then begin
                    colour := thePolys[intersect_polygon-1].polyColor;

                    MakeViewRay(illumPoint, closestIntersectionPt.x3D, closestIntersectionPt.y3D, closestIntersectionPt.z3D, lightRay);

                    for i:= 1 to numPolys do begin
                        if (intersect_polygon <> i) & Intersect(lightRay, thePolys[i-1], my_storage^.poly_extra[i-1], thisDistance, intersectionPt) then begin
                            colour.red := band(colour.red, $0FFFF) div 2;
                            colour.green := band(colour.green, $0FFFF) div 2;
                            colour.blue := band(colour.blue, $0FFFF) div 2;
                            leave;
                        end
                    end;

                    LongintPtr(pixels)^ := bsl( band(colour.red, $0FF00), 8 ) +
                            band(colour.green, $0FF00) + bsr( band(colour.blue, $0FF00), 8 );
                end else begin
                    LongintPtr(pixels)^ := 0;
                end;
                longint(pixels) := longint(pixels) + 4;
            end;
            longint(pixels) := longint(pixels) + rowbytes_add;
        end;
    end;

end.